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ABSTRACT 

Kipp, Carl R. MSME . , Purdue University, August 1985. 
Prediction of Sound Fields in Acoustical Cavities Using 
the Boundary Element Method. Major Professor: R. J. 
Bernhard, School of Mechanical Engineering 

A method is developed to predict sound fields in 
acoustical cavities. The method is based on the indirect 
boundary element method. An isoparametric quadratic boun- 
dary element is incorporated. Using this method, either 
pressure, velocity and/or impedance boundary conditions 
may be applied to a cavity. The capability to include 
acoustic point sources within the cavity is implemented. 

The method is applied to the prediction of sound 
fields in spherical and rectangular cavities. All three 
boundary condition types are verified. Cases having a 
point source within the cavity domain are also studied. 
Numerically-determined cavity pressure distributions and 
responses are presented. The numerical results correlate 
well with available analytical results. 


CHAPTER 1 


INTRODUCTION 


Sound fields in cavities affect many aspects of 
everyday life. It is well known that distribution of 
sound in an interior work environment can markedly influ- 
ence the performance of those working within that environ- 
ment. Most modern modes of transportation place the trav- 
eler in some form of cavity. The perceived comfort level 
of the traveler may again be influenced by the acoustical 
characteristics of the transporting cavity. Consequently, 
noise control engineers in the architectural and transpor- 
tation industries become involved in optimizing acoustical 
cavities so as to minimize the possible negative effects 
on persons residing within those structures. 
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Computer aided design (CAD) procedures are very bene- 
ficial to individuals participating in noise control 
activities. Using CAD procedures, the acoustic charac- 
teristics of a cavity can be determined analytically, usu- 
ally with less time and expense than required for 
experimentally-based techniques. Noise source identifica- 
tion can be accomplished with computerized modeling tech- 
niques. In addition, CAD procedures can readily be app.lied 
in optimizing the acoustical characteristics of a cavity 
after the noise sources have been identified. 

The objective for the research presented here is the 
development of a numerical tool based on boundary element 
theory to predict the sound field within an acoustical 
cavity. The theory incorporates an indirect boundary ele- 
ment method utilizing isoparametric quadratic boundary 
elements. In addition to the usual boundary conditions of 
acoustic pressure and velocity, impedance boundary condi- 
tions are also formulated and implemented. Furthermore, 
capability to model acoustical point sources within the 
cavity domain is included. With all of the aforementioned 
features, the computer program developed during this 
research will provide a very useful tool for noise control 
activities . 

The content of this thesis is organized in the fol- 
lowing manner. Chapter 2 contains an extensive literature 


review on numerical predictions of sound fields. The 
theoretical background of the indirect boundary element 
method for acoustics is presented in chapter 3. Numerical 
implementation of the indirect boundary element theory is 
presented in chapter 4. Chapter 5 presents case studies 
completed through the use of the computer program result- 
ing from this research. Concluding remarks are given in 
chapter 6 . 


CHAPTER 2 


LITERATURE REVIEW 


Many researchers have been investigating techniques 
for numerical prediction of acoustic fields. Consequently, 
a significant body of technical and academic literature 
has been published. In contrast to numerical techniques 
developed for specific problems, such as Green's function 
methods and series solutions, the emphasis in this work 
will be generalized techniques. Such techniques are versa- 
tile and powerful and will someday play an important role 
in computer-aided noise control design. This chapter 
presents a review of the literature applicable to the 
numerical prediction of sound fields in acoustical cavi- 
ties using generalized techniques. Two methods are used 
primarily, the finite element method and the boundary ele- 
ment method. The existing applications of each will be 


discussed. 


2.1 finite Element Methods in Cavity Acoustics 


The finite element method is the most widely-known 
numerical tool available for engineering analysis of com- 
plex structures. A recent survey of finite element tech- 
nology [1] contains a list of 155 textbooks and monographs 
relating to finite element technology. This list does not 
include the numerous journal articles relating to the 
topic. It is obvious from the survey that the finite ele- 
ment method (FEM) is now firmly rooted in various fields 
of engineering analysis. 

In acoustics, the finite element method is an esta- 
blished analytical tool. Most commercially-available fin- 
ite element codes (e.g., NASTRAN, ANSYS, SAP IV) are 
designed to be flexible, thereby providing the ability to 
solve various classes of problems. With an understanding 
of the specific code being utilized and the necessary 
acoustical f undamentals , a finite element code may be used 
to solve acoustical problems as shown by Bernhard [ 2 ]. 

The majority of finite element acoustical studies can 
be classified into two primary categories based on the 
application. These two categories are (1) duct acoustics 
and (2) cavity acoustics for automotive applications. 


2.1.1 Duct Acoustics 


A common analytical method used to analyze acousti- 
cal duct systems is four -pole network theory. In four- 
pole network theory, a transmission matrix can be used to 
define the relationships between the input and output 
variables of acoustic pressure and particle velocity for 
an acoustic device (e.g., a pipe, an expansion chamber , a 
Helmholtz resonator). The sound attenuation or transmis- 
sion loss of a device may be determined using the 
transmission matrix. The transmission characteristics of 
an acoustical system comprised of a series of well-defined 
devices can be found by combining the respective transmis- 
sion matrices. 

Young and Crocker [3] employed the finite element 
method to determine the four-pole parameters of transmis- 
sion matrices for expansion chambers. Since the initial 
studies were for simple expansion chambers, two- 
dimensional, four-node rectangular elements were used. 
Each element had twelve degrees-of-f reedom : a pressure 
and two pressure gradients, one in each coordinate direc- 
tion, at each of the nodes. Utilizing fourth-order Hermi- 
tian shape functions, the transmission loss of a rigid- 
wall expansion chamber was numerically evaluated and found 
to correlate well with theoretical results. The chamber 
was driven at the inlet with a harmonic velocity and 
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terminated at the outlet first with a rigid-wall boundary 
condition and then with a pressure release boundary condi- 
tion. 


As an extension of this work, Young and Crocker stu- 
died complicated muffler systems having multiple chambers 
[4]. Furthermore, the effects of non-rigid walls were 
included. Numerically, a multiple chamber system is 
assimilated by the cascade connection of the transmission 
matrices for the individual components. For a more com- 
plex system, the finite element method was well-suited for 
determining the four-pole parameters of the individual 
components to be used in calculating the transmission loss 
for the total system. 

Astley and Eversman examined acoustical propagation 
abcve the cut-off frequencies in ducts in terms of the 
eigenvalue problem [5]. The research included the effects 
of steady mean flow or sheared flow within the duct. Capa- 
bility was incorporated to study lined ducts by consider- 
ing wall admittance. Two types of elements were con- 
sidered. The primary element was a three-node element 
with quadratic shape functions. The variables in the for- 
mulation were the perturbed velocities in the x and y 
directions and the perturbed pressure. With the primary 
element, extraneous eigenvalues were found to result at 
relatively high frequencies. It was felt that the 
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extraneous eigenvalues were due to slope discontinuities 
in the quadratic shape functions. To overcome this diffi- 
culty, two-node elements incorpor acing shape functions 
having slope continuity were considered. The extraneous 
eigenvalues were eliminated although the solution accuracy 
was reduced at the higher eigenvalues. The finite element 
method proved to be applicable to the solution of an 
eigenvalue problem for ducts. 

In addition to using acoustical finite elements, 
Craggs studied the use of absorptive finite elements for 
acoustical applications [6]. The implicit advantage of 
absorptive elements is that extended reaction of linings 
can be included. Most techniques supply the lining admit- 
tance as a boundary condition thereby only allowing for 
local reaction of the lining. Using the absorptive ele- 
ments, the absorption coefficients for various resistivi- 
ties of a constant thickness material were determined and 
found to compare quite well with analytical results. The 
acoustic and absorption finite elements were also coupled 
to study a lined expansion chamber. The transmission loss 
of the chamber as a function of lining resistivity was 
determined. Results for the lined expansion chamber would 
be very difficult to determine analytically. Thus, the 
advantage of the absorption finite element is evident. 
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The finite element method may also be used to examine 
the sound field within a duct as was shown by Cederfeldt 
[7]. Towards that objective, Cederfeldt studied a right- 
angle bend, a lined rectangular duct and an expansion 
chamber. The models were two-dimensional and comprised of 
four -node quadrilateral acoustic elements. Contour maps 
of the sound fields within the structures were presented. 
A variation of the right-angle bend created by adding a 
lining on the exterior corner provided quite a different 
sound field as compared to the bend without the lining. 
Even with the somewhat simplistic elements, dramatic 
representations of the sound fields were obtained. 

A more sophisticated element was used by Ling, Hamil- 
ton and Allen to determine the sound field in axisymmetr ic 
ducts [8]. The element was a two-dimensional, iso- 
parametric, axisymmetr ic element with cubic polynomials as 
the shape functions. With this element, the sound field in 
a bottle-lifce duct was determined numerically and compared 
with one-dimensional theoretical results and experimental 
measurements. Two types of lined ducts were also studied: 

(1) a duct with segmented lining (i.e., hard wall and 
impedance boundary conditions on different segments) and 

(2) a duct with an exponential distribution of absorbent 
material along its length. Both of the above lining condi- 
tions produce a pressure distribution within the duct 
which is difficult to obtain analytically but, as shown by 
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these researchers, can be obtained with the finite element 
method . 

Whereas Ling, Hamilton and Allen used an iso- 
parametric element, Doyle and Faulkner applied a sub- 
parametric element for acoustic pressure analysis [9]. The 
geometry variation within the element was linear while the 
pressure variation was made to be cubic under the philoso- 
phy that the subparametr ic element would be more computa- 
tionally efficient than an isoparametric element. At each 
grid point, the four independent variables of pressure and 
the pressure gradient in each of the three orthogonal 
coordinate directions were used. An additional feature 
highlighted in this research is the ability to assign dif- 
ferent fluid densities and speeds of sound at any grid 
point. With the method as described above, the eigenvalue 
problem for a uniform hard-walled duct was solved to find 
both the natural frequencies and the mode shapes of the 
duct. A second type of problem was analyzed where a har- 
monic velocity was assigned to one wall of the duct. The 
pressures in the duct were determined for a hard-walled 
case and cases with acoustic impedances on one or more 
walls . 

2.1.2 Cavity Acoustics for Automotive Applicatio ns 

Cavity acoustics problems for automotive applica- 
tions are generally more complex than those for duct 
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acoustics problems. The geometry of the structures are not 
simple rectangular or cylindrical shapes. In addition, the 
interior of the cavity may contain obstacles affecting the 
sound field. Forcing functions can not normally be defined 
with ease when forced response problems are considered. 
For these reasons, the research presented below has been 
separated from the duct acoustics research presented in 
the previous section. 

Craggs studied the sound field in a passenger car 
cavity when the cavity is coupled to the engine compart- 
ment through a flexible boundary [10]. The boundary 
motion was formulated in terms of the normal modes of the 
enclosure. The governing equations for a single enclosure 
with a flexible boundary and the structural coupling 
between two enclosures were developed. An example of a 
passenger car compartment coupled through a flexible boun- 
dary with the engine enclosure was then studied using the 
numerical formulation. The model could be used to study 
the sound field induced within the cavity due to an exci- 
tation within the ' engine compartment. From the results, 
some generalizations were made regarding the situations in 
which the passenger cavity either was or was not greatly 
influenced by the excitation in the engine compartment. 

Sung studied finite element applications to the 
acoustics of a passenger car compartment and an engine 
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combustion chamber [11]. The NASTRAN finite element code 
was used for both cases. Three-dimensional models of the 
passenger compartment were constructed for two cases: (1) 
a compartment with a bench seat and (2) a compartment 
with bucket seats. The resonant frequencies of the cavi- 
ties were determined and the corresponding mode shapes 
presented with an acoustic isobar representation. Inclu- 
sion of the seat in the model provided a more realistic 
situation than a model without the seat. Por the combus- 
tion chamber model, a piston head with an asymmetric bowl 
was included in the chamber model. Again, the resonant 
cavity frequencies and mode shapes were determined from 
the finite element model. These resonant frequencies were 
compared with experimentally determined resonant frequen- 
cies and found to have good agreement. The final observa- 
tion was that the acoustic finite element models were suc- 
cessful in predicting cavity boom frequencies in passenger 
compartments and knock- induced frequencies in a combustion 
chamber . 

In conjunction with their work in duct acoustics, 
Doyle and Faulkner extended the research by studying the 
sound field in an automobile interior [9]. The first ana- 
lyses performed were similar to those done by Sung [11] in 
that the resonant frequencies and mode shapes were quanti- 
fied. However, more interesting results were obtained for 
a forced response problem. Experimentally determined 


pressures at * sun roof opening were supplied to the fin- 
ite element model as the forcing function. The acoustic 
field produced by the forced model was compared with a 
discrete number of points at which experimental data had 
been collected. Although the numerical pressures were 
slightly less than the experimental pressures,, there was 
good correspondence between the two sets of data. 

The above-cited references exemplify the types of 
acoustic cavity problems which may be studied with finits 
element techniques. Moreover, they indicate the popular- 
ity of finite element acoustics. The following section 
presents acoustical applications of the Boundary Element 
Method. Although not yet as popular as finite element 
techniques, boundary element techniques share many of the 
attractive features of the finite element method and are 
probably better adapted to certain problems. 

2.2 Boundary Element Acoustics 

There are two basic classifications of Boundary Ele- 
ment Methods: (1) Direct Boundary Element Methods (DBEM) 
and (2) Indirect Boundary Element Methods ( IBEM) . The dis- 
tinction between the two methods must be made at this 
point for clarification of the following discussion. Por 
the current project, an IBEM formulation is utilized as 
will be shown in Chapter 3. 
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The difference between the two class if icati*-ns is 
evident in the integral equations. For the DBEM, "the unk- 
nown functions appearing in the integral equations are the 
actual physical variables of the problem' [12]. The physi- 
cal variables for an acoustics formulation are the acous- 
tic pressure and the particle velocity. A general form of 
the DBEM integral equation is 

p(€) - jTp (x,£)v(x) - p ( x ) v (x,*)ldB 
B 

+ /p*(x, T)F , (x)dD (2.1) 

D 

where p(4) is the pressure at any point, p(x) is the pres- 
sure at a boundary point, v(x) is the normal velocity at a 

* — — * — — 

boundary point, p (x,£) and v (x,4) are the respective 
fundamental preesure and velocity functione from a boun- 
dary point to any point and ?(x) represents applied source 
strengths over the domain of the problem. 


The formulation for an IBEM is significantly dif- 
ferent from that for a DBEM. For the IBEM, 

"The integral equations are expressed entirely in 
terms of a unit singular solution of the original 
differential equations distributed over the boun- 
daries of the region of interest. The density func- 
tions themselves have no specific physical signifi- 
cance tut once they have been obtained from a numeri- 
cal solution of the integral equations the values of 
the solution parameters anywhere within the body can 
be calculated from them by simple integration 
processes. " [ 12] 

Or, in other words, the boundary is replaced by a distri- 
bution of sources which reproduce the specified boundary 
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solution and allow computation of the physical variables 
anywhere in the problem domain or boundary. Separate 
equations are required for portions of the boundary 
depending upon the specified boundary conditions on the 
particular portion of the boundary. The general form of 
the IBEM integral equation for a pressure boundary condi- 
tion is 

p(7) - /o(x)p*(x,7)dB + /p*(x,7)*(x)dD (2.2a) 

B D 

whereas the general equation for a velocity boundary con- 
dition is 

V(7) - ±co(x) + /o(x)v* (x,7)dB + Jv*(x,7mx)dD(2.2b) 

B D 

where o(x) represents the fictitious source density func- 
tions at the boundary points and c represents an integra- 
tion constant due to integral singularities which will be 
further explained in chapter 3. 

Although the details of the two formulations are not 
critical at this point, an understanding of the existence 
of differing formulations is important. Some of the fol- _ 
lowing reviewed works use a DEEM while others use an IBEM. 
This point will be emphasized within each review. 
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2.2.1 I nyggt ia at ifl ng Qf Exterior Acoustical ProhUma 

Most of the boundary element research heretofore 
published has dealt with the acoustic field of a domain 
exterior to a vibrating or scattering structure. Pri- 
marily, this is a consequence of the capability of boun- 
dary element methods to model domains extending to infin- 
ity, a task not easily accomplished with finite element 
models. However, the integral formulations of interior 
regions are nearly equivalent to those for exterior 
regions. Hence, much of the development for exterior 
regions may be extrapolated to developments for interior 
regions . 

Chen and Schweikert [13] were the earliest contem- 
porary researchers to apply boundary element techniques to 
acoustic radiation problems. A numerical implementation of 
Huygen’s principle, which is essentially a primal IBEM 

V 

formulation, was utilized whereby a distribution of simple 

* 

sources was imagined over the bounding surface. The sur- 
face itself was considered to be an array of triangular 
elements over which a uniform velocity with a constant 
phase angle was prescribed on each element. Each element 
was analogous to a rigid piston vibrating with a constant 
harmonic behavior. After defining a model by its geometry 
and boundary conditions, the boundary integral equations 
were numerically evaluated to calculate the simple source 
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distribution over the model boundary. The field pressures 
were quantified from the source distribution on the boun- 
dary. Two examples were given. The first example consisted 
of a piston in a rigid sphere for which the exterior field 
pressure radiation pattern was calculated. Likewise, simi- 
lar exterior pressure patterns were determined for a stif- 
fened cylinder immersed in water. The fluid-structure 
interface effects were included for the second example. 

A DBEM was incorporated into a study of acoustic 
radiation from surfaces of revolution by Chertock [14]. 
The variables in the integral equations for a DBEM 
represent the acoustic pressure and particle velocity. So 
by specifying one type of boundary condition, the remain- 
ing unknown boundary condition may be determined with the 
Helmholtz integral equations. Furthermore, the field pres- 
sures can be calculated based on the boundary parameters. 
Chertock used this procedure to numerically evaluate the 
surface pressures and field pressures from the velocity 

boundary conditions for rigid-body vibration of a sphere, 

# 

quadrapole vibration of a prolate spheroid (i.e., a 
cigar-shaped surface) and quadrapole vibration of an 
axisymmetric surface which was unsymmetric about all axes 
except tiie major axis. Surfaces of revolution were stu- 
died because theoretical results may be obtained for com- 
parison with numerical results. Good correlation was found 
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between the numerical and theoretical solutions for the 
cases in question. 

Copley [15] also applied a DBEM to the study of 
acoustic radiation from axisymmetric radiators. However , 
Copley [15] used an approach which differed from that of 
Chertock [14]. Copley utilized the Interior Helmholtz 
integral equation whereas Chertock utilized the Surface 
Helmholtz integral equation. In the interior Helmholtz 
integral approach, the field points of integration lie 
within the boundary of the vibrating body and, in fact, 
are located along the axis of symmetry for axismmetric 
radiators. The corresponding field points lie on the 
boundary in the surface Helmholtz integral approach. 
Hence, the appropriate titles of each approach refer to 
the location of the field integration points. Using the 
Interior Helmholtz integral equation, the acoustic radia- 
tion from a spheroid and a finite cylinder with a capped 
end were studied. The farfield directivity pattern of the 
finite cylinder was compared to results obtained by Willi- 
ams, et.al. using a method based on expansion in spheri- 
cal harmonics [16]. Generally, the results correlated 
well with the exception of some discrepancy which was 
thought to be due to differences in which the velocity 
distributions were prescribed. 


Taking a more general point of view, Copley [17] also 
studied failures which oc.cur with integral formulations 
for acoustic radiation problems. At certain frequencies, 
the integral formulations are unable to represent the 
steady-state harmonic radiation from a finite, smooth, 
closed surface on which normal velocity is prescribed. 
These failures occur at the interior Dirichlet eigenfre- 
quencies for both the IBEM and the Surface Helmholtz 
integral equation which is a DBEM. Copley [17] presented 
proofs that the failures were due to nonexistence of 
sources for the IBEM and as nonuniqueness of sources and 
doublets for the Surface Helmholtz integral equation. The 
failures are inherent in the formulations and do not 
represent physical attributes of the problem. However, 
the Interior Helmholtz integral formulation used by Copley 
[15] does not suffer from these failures. 

Proofs similar to those of Copley [17] were presented 
by Burton and Miller [18]. Moreover, suggestions to over- 
come the failures were made. One alternative consisted of 
over determining the system of equations by using both the 
Surface Helmholtz integral equations and its differen- 
tiated form. The second but more economical alternative 
was to combine the two types of equations into a single 
set of equations which would not overdetermine the system. 
No computational examples were provided for these 
approaches . 
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A combined Helmholtz integral equation formulation 
was developed by Schenck [19] to overcome the deficiencies 
of the above mentioned formulations. The technique ini- 
tially applyed the Surface Helmholtz integral equations 
and then overdetermined the system of equations with the 
Interior Helmholtz integral equations. The overdetermined 
set of equations was solved by a least squares orthonor- 
malizing procedure to evaluate the surface pressures. From 
the surface pressures, the field pressures were calculated 
with the Helmholtz integral equation. Essentially, the 
described procedure is a DBEM as all the equations are of 
the direct formulation type. A number of numerical exam- 
ples were presented for the combined Helmholtz integral 
equation formulation of which only three will be dis- 
cussed. The first example is a uniformly vibrating sphere. 
The problem illustrates the advantage of the combined 
equation formulation. At a characteristic wave number, the 
Surface Helmholtz integral equations fail to yield the 
correct surface pressure. However, the combined Helmholtz 
integral formulation was capable of correctly solving for 
the surface pressure on the sphere by specifying only one 
point within the sphere. A second example of a right cir- 
cular cylinder with rigid ends on which velocity boundary 
conditions were prescribed was studied. The resulting far- 
field pressure pattern was found to be consistent with 
that obtained by Chertock [14] but inconsistent with that 
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obtained by Williams, et.al. [16]. As a variation on this 
model, the velocity distribution was revised to reflect 
the velocities that would occur at the boundary if a 
source were at the center of the cylinder and the cylinder 
was not physically present. The acoustic pressure radia- 
tion pattern was equivalent to a pattern for a simple 
source as was expected. A rectangular parallelepiped was 
the third model considered. Again, the velocity distribu- 
tion due to a simple source was assigned on the boundary. 
The resulting radiation pattern was the expected simple 
source radiation pattern. For both of the latter cases, 
the Surface Helmholtz integral equation without the inte- 
rior points failed at the characteristic wave numbers. 

The concept of overdetermination of the system equa- 
tions was also addressed by Paiszczyk and Klosner [20]. 

l 

The overdetermining equations were developed with the 
Exterior Helmholtz integral equations unlike Schenck's 
research which used the Interior Helmholtz integral equa- 
tions for the overdetermination [19]. An iterative pro- 
' cedure is required in the method proposed by Piaszczyk and 
Klosner [20]. At the surface, an approximate acoustic 
fluid impedance was assumed. The Exterior Helmholtz 
integral equations were used to solve for approximate 
pressures at selected field points. The approximate field 
pressures were then applied in conjunction with the Exte- 
rior and Surface Helmholtz integrals to solve for the 
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surface pressures. A least squares procedure was required 
to calculate the surface pressures. The surface pressures 
were then substituted into the overdetermined set of equa- 
tions to solve for the field pressures which were con- 
sidered to be the next iterative value of the field pres- 
sures. The entire iterative process continued until a con- 
vergence was established. For this exterior overdetermina- 
tion scheme, the selection of the overdetermining points 
is not critical due to its iterative nature. Thus, the 
scheme may be applied without loss of generality to struc- 
tures of arbitrary shape. A number of computational exam- 
ples were presented displaying the versatility and accu- 
racy of the exterior overdetermination method. 

The research of Schenck [19] was used in combination 
with finite element analysis to study acoustic radiation 
from sonar transducers by Smith, Hunt and Barach [21]. 
Although the combined Helmholtz integral equation formula- 
tion was applied, only the Surface Helmholtz integral 
equations were utilized in the sonar transducer study. A 
specific acoustic impedance matrix was derived with the 
integral equations. Definition of consistent mass and 
stiffness matrices for the complex structures came from 
the finite element analysis. The total structural- 
acoustical system of equations combined the mass and 
stiffness matrices with the acoustic loading included as 
surface loading forces. Finally, the pressure radiation 
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patterns were determined using the results of the 
structural-acoustical analysis and the Helmholtz integral 
equations to evaluate the pressures at points in the field 
surrounding the transducer. Numerical results for a 
piezoelectric sphere and a piezoelectric free-flooded 
cylinder were compared with experimental results and found 
to agree within 5 per cent. Some of the error was attri- 
buted to the piezoelectric nature of the material in the 
transducers . 

Engblom and Nelson [22] expanded upon the work of 
Smith, Hunt and Barach [21] by the addition of two 
features. The first additional feature involved allowing 
quadratic variation of the acoustic variables and linear 
variation of the geometrical variables over the surface of 
an element. Prior to this work- a constant variation of 
the acoustic variables had been assumed. In doing sc f 
discontinuities of acoustic variables may arise at the 
boundary of neighboring elements. A quadratic variation 
ensures that the acoustic variables become continuous at 
the element boundaries. The second feature added was a 
coordinate transf ormation used to integrate around the 
singularity which exists in the Helmholtz integral formu- 
lations. Since the elements were triangular, the coordi- 
nate transformation involved a change from area coordi- 
nates to polar coordinates. In doing so, the singularity 
was removed thus allowing the integration to be completed 



i , without the singularity complications. The common examples 

of acoustic radiation from a sphere and a right circular 
cylinder were numerically evaluated with good results. 

Unlike the approach originally developed by Schenck 

\ 

[19] to overcome the difficulty of nonuniqueness, Meyer, 
Bell, Zinn and Stallybrass [23] implemented an approach 
suggested by Burton and Miller [18]. That is, the system 
of equations consisted of the Surface Helmholtz integral 
equation and its differentiated form, both of which were 
direct formulations. The system of equations was not an 
overdetermined set but a combination of the two types of 
integral equations. A difficulty arose in that the dif- 
ferentiated form of the Surface Helmholtz integral equa- 
tion contained a strong singularity which could not be 
directly numerically integrated. Nevertheless, the authors 
proved through riiorous mathematics and computational con- 
siderations that the integration was possible. For an 
example of a piston set in a rigid sphere, the error for 
the numerical farfield pressures remained less than 10 per 
cent as compared with exact analytical solutions. 

One further refinement of the DBEM for acoustic radi- 
ation problems was achieved by Seybert, Soenarko, Rizzo 
and Shippy [24]. The authors provided further sophistica- 
tion in the discretization process through the use of an 
isoparametric element. With the isoparametr ic element, 

i 

i 

t 

I 

} 
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both the geometry end the ecouetic variables were interpo- 
lated with quadratic shape functions. Although not immedi- 
ately apparent, the isoparametric element is somewhat more 
refined then the element used by Engblom and Nelson [ 22 ] 
which ussd a linear geometrical variatirn. Consequently, 
efficiency and accuracy were increased as less elements 
were required to model curved or other irregular surfaces. 
To illustrate the isoparametric element capability, the 
examples of a pulsating sphere and an oscillating sphere 
were considered and found to compare well with theoretical 
results except at the characteristic wavenumbers. The 
nonuniqueness of the formulation was not accounted for by 
these authors. 

The preceding reviews in this section represent the 
advances which have occurred in boundary element technol- 
ogy for exterior domain acoustics, particularly for radia- 
tion problems. Many of the same authors have also studied 
scattering problems. Thus far, the theoretical implica- 
tions have been discussed. The remaining portion of this 
section will review applications of the aforementioned 
theory to somewhat more complicated problems. 

A DBEM has been used by Se2nec [25] to study the dif- 
fraction of sound around barriers. The barriers were any 
general type or shape used for noise abatement (e.g., a 
barrier around a roadway to hinder traffic noise). The 
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capability to include reflectivity and absorptivity 
effects of both the ground and the barrier was achieved 
through the application of impedance boundary conditions 
on the surfaces. Prior to the work of Seznec, researchers 
had been unable to adequately account for the reflectivity 
and absorptivity effects. The diffraction of sound around 
barriers of various shapes were determined. For one bar- 
rier, that of a infinitely thin reflecting barrier, the 
numerical solution was favorably compared with solutions 
already existent in the literature. Overall, the DBEM was 
found to be effective as a design tool in studying the 
diffraction of sound around barriers. 

Sadek and his coworkers developed a procedure to 
predict the acoustic emission from a machine using the 
design drawings for the machine [26]. A finite element 
model of the machine structure was used to predict the 
structural modes of vibration and the modal response to 
typical operating loads for the machine. The modal 
response yielded the surface velocities which were used as 
boundary conditions for a DBEM formulation to predict the 
acoustic emission. Two applications of the procedure were 
presented. The first example was for a forging machine 
structure under impulse loading characteristic of the 
machining operation. Reflections of sound from the ground 
could be and were included in determining the acoustic 
emission. Comparison of the numerical results with 
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experimental results shewed that each exhibited the same 
general trends and magnitudes although there were varia- 
tions at specific points. second application involved 
studying the acoustic emission from one component of a 
hydraulic hammer [27]. The component considered to be the 
most offending noise source was studied. The acoustic 
emission analysis procedure was applied to the initial 
design and possible redesign configurations. As a result 
of these analyses, one redes icin configuration quantita- 
tively produced less acoustic emission than the initial 
design and the other redesigns. The process followed in 
this work emphasizes the manner in which a boundary ele- 
ment method may be usod in the design process. 

A method to predict acoustic intensity on the surface 
of a vibrating body utilizing a DBErf formulation has been 
developed by Benner and Koopmann [28,29], Much like the 
previously presented research, a finite element code was 
used to define the dynamics of the structure. Since deter- 
mination of acoustic intensity requires both the pressure 
and velocity, both of these quantities were calculated 
with the boundary element computer code. In addition, 
acoustic pressure and velocity, and hence acoustic inten- 
sity, were calculated at any field point. The boundary 
element acoustic intensity method was applied to the study 
of crack propagation along the boundary of a flat plate 
[30]. To simulate crack propagation, successive collinear 
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degrees -of -freedom were released within the finite element 
model. The acoustic surface intensity was calculated for 
various modes of vibration as each additional degree-of- 
freedom was released. A definite change in the intensity 
distribution occurred as the crack propagated. This appli- 
cation highlighted the possible use of a boundary element 
program as a tool in analyzing cracked or flawed struc- 
tures for non-destructive testing. A second application 
was to study the radiation characteristics of a slender 
box-type structure representative of a box girder in a 
rapid transit transportation system [31]. The intent was 
to isolate high radiation efficiency modes so that damping 
treatments could concentrate on reducing those modes. 
Relatively high radiation efficiency modes were identified 
using the boundary element formulation but analysis of the 
structure with damping was not presented. This application 
highlights one manner in which boundary element formula- 
tions- may be used for design modifications. 

2.2.2 Investigations of Interior Acoustical Problems 

The preceding section and the vast majority of 
literature relating to boundary element technology for 
acoustical applications consider exterior domain problems. 
The following discussion reviews literature applied to 
acoustical problems in the interior domain. 
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Tai and Shaw [32] have used a DBEM formulation to 
determine eigenvalues and eigenmodes of the Helmholtz 
equation for arbitrary domains. The technique is applica- 
ble to both two- and three-dimensional domains. However, 
the numerical example presented was for a two-dimensional 
domain. Within an element, the acoustic variable was con- 
sidered to be constant. The numerical example presented 
consisted of a right triangle along whose hypotenuse the 
eigenmodes were determined after the eigenvalues had been 
found. The first two modes were found for three cases hav- 
ing varying magnitudes of the acute angles. This problem 
was jo chosen to emphasize that the integral formulation 
is capable of providing results when a separation of vari- 
ables approach is not. The separation of variables tech- 
nique can be used when the boundaries are along constant 
coordinate lines thereby allowing the equations to 
separate into individual coordinate directions. For the 
case of a right isoceles triangle, which has a known 
analytical solution, the numerical solution with the DBEM 
formulation and analytical solution showed excellent 
agreement. 

Tanaka, Fujikawa, Abe and Utsuno developed a method 
incorporating a DBEM formulation to calculate the transfer 
matrices used in analyzing muffler systems [33]. This 
method was found to yield good accuracy and better effi- 
ciency than the finite element approach applied by Young 


and Crocker [3,4]. Boundary element theory was applied in 
regions where it is difficult to apply plane-wave theory. 
Plane-wave theory was applied in regions where that theory 
was still applicable. The transfer matrix for the entire 
system was obtained by multiplying the matrices from each 
component and used to calculate the insertion loss. 
Impedance boundary conditions were included for some 
cases. The method combining boundary element and plane- 
wave theory gave good results and also decreased both 
modeling and computational time as compared to other 
numerical techniques. 

Sestieri, Del Vescovo and Lucibe.Ho have considered 
the problem of structural-acoustical coupling in cavities 
by application of a DBEM [34]. The effects of the acoustic 
loading were included as forces on the structure in an 
inhomogeneous Helmholtz formulation (i.e., a forced 
acoustical problem). An overdetermination procedure simi- 
lar to that used by Schenck [19] was used whereby the 
resulting pressures were a least squares solution of the 
problem. Constant acoustic variable variation was implied 
on an element surface. Impedance boundary conditions were 
also included. A series of studies involving a cylindrical 
cavity excited by a harmonic driver at one end produced 
very good agreement with experimentally obtained results. 
Other studies investigating the quantity and placement of 


the interior overdetermining points were also presented 
but will not be discussed here. 


2.3 Conclusions 

This chapter has presented a number of references 
related to prediction of sound fields in cavities and 
boundary element applications in acoustics. It reviews the 
status of techniques being applied to the problems within 
the concern of the current research. In addition to the 
technical literature, there are good academic references 
available relevant to boundary element techniques for both 
general and specific topics [12,35,36]. 

A few developments from the literature review should 
be emphasized due to their significance to the current 
research. IBEM's require less information than DBEM's. 
Equations 2.1 and 2.2 illustrate that less information is 
required in the integral equations. More importantly, only 
one quantity, the surface source distribution, is ini- 
tially computed at the boundary whereas, for DBEM's, both 
pressure and its gradient, depending on the boundary con- 
dition, are resultants. Hence, the IBEM may be prefered 
for noise source identification since interpretation of 
the results would be more straightforward. Impedance 
boundary conditions have been applied in past research and 
will be applied in the current research. The information 
relating to the failures of the boundary element 
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formulations may provide guidance in case difficulties of 
this nature are encountered. The literature review has 
shown that isoparametr ic representation, which is to be 
utilized in the current research, has been previously 
utilized successfully. One very important concept 
discovered during this review is that of domain transfor- 
mation to integrate around the singularity inherent in the 
Helmholtz integral formulations. A similar approach will 
be used and is detailed in chapter 4. Finally, the litera- 
ture review shows a definite lack of boundary element 
applications to acoustical cavity problems. Of the 21 
boundary element references cited, only 3 were concerned 
with cavities. 


CHAPTER 3 


THEORY 


The problems to be considered are those dealing with 
acoustical behavior within cavities. A description of the 
problem is shown in Figure 3.1. A point within the domain 
D is located by the vector x. Likewise, a point on the 
boundary B is defined by the vt ^tor <. The vector normal 
to the boundary is signified by n and is considered to be 
directed outward from the cavity. An applied distributed 
source of strength t within the domain is located by the 


vector x . 

3 



Por a problem of this type, the governing equation is 
the familiar non-homogeneous Helmholtz equation. 

2 2 — 

v • + k • - f ( Xg ) (3.1) 

where x fl is the vector locating the applied source. Equa- 
tion 3.1 is the linearized, lossless Helmholtz equation 
formulated in terms of velocity potential, <*>. The vari- 
ables of acoustic pressure and particle velocity can be 
related to the velocity potential through the relation- 
ships [37] 


p - ~jwp4> (3.2) 
v - TO. (3.3) 

A Huygen's principle formulation which may be used to 
solve the homogeneous form of equation 3.1 has been 
developed by Chen and Schweikert [13]. In this formula- 
tion, a distribution of simple sources is considered to 
exist at the boundary of the cavity in question. The boun- 
dary source magnitudes are determined in such a way that 
the boundary conditions for the problem are satisified. 
Chen and Schweikert' s method is essentially a simple 
indirect boundary element method. ^ o develop a more com- 
plete indirect boundary element method, a variable distri- 
bution of sources may be assumed as will be shown in 
Chapter 4. The velocity potential at any point due to the 
assumed source distribution is 
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e -jkr 

®(x) - Jo K C )~ — dB (3.4) 

B 

where r is the distance between the points x and C and 
0(<) is the source strength. Although equation 3.4 is 
written for a point in the domain, it applies equally as 
well for a point on the boundary. Similarly, the gradient 
of the velocity potential is 

TO(x) - 4* c o(x) + S °(O v 

3~ 

where c results from a singularity integration which will 
be defined further in this chapter and B~ is the boundary 
excluding a small region at the singularity. 

At this stage, a definition of two quantities will be 
made to ease notational difficulties. Combination of equa- 
tions 3.2 and 3.4 allows an integral representation of the 
pressure at a point to be expressed as 

e~ jkr 

P(x) - -jwpjocc) — - — dB. (3.6) 

B 

* 

A quantity known as the fundamental pressure solution, p , 
is defined by 

t -jkr 

p ( C , x ) - -jup 5 -^ . (3.7) 

Hence, equation 3.6 may be more concisely written as 

PU) - Jo(()t)*U,x) dB. 


,-jkr 


dB 


(3.5) 


(3.8) 
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* 

Likewise, the fundamental velocity solution, v , may 
be defined by considering, the partial derivative within 
the integral portion of equation 3.5, that is 




(3.9) 


Since, for a well-posed boundary value problem, the com- 
ponent of velocity of interest is the normal velocity, the 
following relationship from Kaplan [38] is used. 


^ - *() ' h (3.10) 

an 

Using equations 3.9 and 3.10, the fundamental velocity 
solution may be defined as 


* — — 


A. 


v (C ,x) - n 'Vr 


zl 

r 2 r 


(3.11) 


Equation 3.5 may now be written more concisely as 


* — — 


v(x) - 4vco(x) + J o(C)v ( C , x) dfl 


(3.12) 


B 

where v is the normal velocity component of the velocity 
potential gradient in equation 3.3. 


The quantity c in equation 3.12 results from integra- 
tion around a singularity point. Prom examination of the 
fundamental solutions, it is evident that as the distance 
between the two points of concern approaches zero the fun- 
damental solutions become infinite. However, the fundamen- 
tal pressure solution is a ^ singularity while the funda- 
mental velocity solution is a — j- singularity. The 
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singularity is weak and may be integrated but the “4 

r z 

singularity is strong and must be treated as a Cauchy 
principal value [40]. The basic premise is to exclude a 
small region surrounding the singularity from the integra- 
tion and account for the exclusion with the added 'free 
term 1 (i.e., the quantity c) [41,42]. Evaluation of the 
’free term' for an interior domain can be done with the 
following rules. 


c - 0 for an exterior point 


c - -1 for an interior point 


c - ~~2 for a point on a smooth boundary 
c - for a point on a corner 


A smooth boundary is any boundary with a unique tangent 
plane. For a corner point, n is the value of the salid 
angle included in the domain at the corner expressed in 
steradians . 

The previous discussion is for the homogeneous prob- 
lem but the inhomogeneou- problem can be solved by aug- 
menting equations 3.8 and 3.12 with a domain integral 
[39]. The domain integral represents the contribution of 
an applied distributed source of strength ?. After inclur 
ina the domain integral, the equations become 
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and 


* — — 




p(x) ■ /o(C)p (C ,x) dB ♦ JV(x«) p (x«,x) dD (3.13) 

B D 8 8 


* — — 


v ( x ) - 47T C a(x) + / o ( C ) v ( < , x) dB + 


B 


JV(x )v (x ,x) dD 

D 8 8 


(3.14) 

where x indicates the location of the applied source. 

s 


If the applied sources are considered to be monopole 
sources, evaluation of the domain integral is signifi- 
cantly simplified. For a point source, the domain integral 
is easy to evaluate as the point source may be represented 
by a Dirac delta function multiplied by the source 
strength. Thus, the integral exists only at the point of 
application of the source. The domain integral becomes 






_ _ nsource _ _ 

JVU B )P (*.'*) dD “ L *f*i a ( x -i ~ *)P <x a ,x) dD 

D 8 8 i-i D 1 81 8 


nsource t _ _ 

L 4^ p (x ai ,x) 
i-1 1 81 

for a pressure boundary condition and 


(3.15) 


_ * _ _ nsource _ _ * _ _ 

St(x )v (x ,x) dD - L /♦ i 0(x_, - x)v (x_,x) dD 
D 8 8 i-1 D 1 81 8 


nsource t _ _ 

L ATrt.y (x 8i ,x) 

i-1 


(3.16) 


for a velocity boundary condition where nsource is the 
total number of applied sources. Distributed internal 


sources are not difficult to evaluate but are rare. 


Vibrating sources which do not behave as point sources can 
be treated as additional boundaries. 


Impedance boundary conditions may also be applied 
with this IBEM. Specific acoustic impedance is defined as 
the ratio between pressure and particle velocity [43] or 

2 - (3-17) 
Manipulation of equation 3.17 yields an equation which is 
applicable to the integral equations. 


p - zv - 0 (3.18) 
In terms of the integral equations, impedance boundary 
conditions may be formulated as 


-4vzc0(x) + /«*(£) [ 

3 

nscurce 

L 4t r* (x ) 

i-1 


P (</*) 


p'<x si .iO 


1 s 

Equation 3.19 can be solvec. for 
source distribution o. 


- zv* (C ,x)j dB + 

zv*(x 8 .,x)] - 0 . (3.19) 
the unknown fictitious 


All the integral equations have now been developed 
for this indirect boundary element formulation. Depending 
on the type of boundary condition applied, equations 3.13, 
3.14 and/or 3.19 are used to solve for the assumed source 
density distribution, o, on the boundary of the body in 
question. Once the source density distribution has been 
evaluated, the same three equations may be applied to 
solve for both the unknown boundary conditions and the 
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acoustic variables at desired points within the domain. 
The numerical implementation of the proceeding theory is 
presented in the next chapter. For the purpose of th* 
current research, only the acoustic pressures at the boun- 
dary and field points will be evaluated as those quanti- 
ties adequately describe the sound field within the cav- 


ity. 


CHAPTER 4 


NUMERICAL IMPLEMENTATION 

The preceding chapter describes the theory of the indirect 
boundary element formulation of the cavity acoustics prob- 
lem. To implement the theory into a usable computer pro- 
cedure, a number of numerical analysis techniques were 
applied. This chapter presents the details of the numeri- 
cal analysis techniques which were used for the implemen- 
tation. 


4.1 Element Definition 


The boundary element incorporated into this research 
is an isoparametric element with quadratic shape func- 
tions. The term isoparametric means the same interpolation 
function is employed for both the geometric approximation 
and the acoustic variable approximation. For an iso- 
parametric element, the number of nodes is equal in both 


the local and the global coordinate systems. Consequently, 
a node in the local coordinate system correlates directly 
to a node in the global system. The local element is shown 
in figure 4.1. 4n example of a global element is shown in 
figure 4.2. 



1 5 2 


Figure 4.’ • Local Isoparametric Element 
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Figure 4.2 - Global Isoparametric Element 


The important feature of the element is the shape func- 
tions. The shape functions are used to interpolate the 
variation of the variables within an element. For the iso- 
parametric element, one set of functions describe the 
variation of both the element geometry and the solution 
variaoles. The element in figure 4.1 has a 3x3 nodal dis- 
tribution in a two-dimensional space (i.e., the $ and the 
V directions). Consequently, a quadratic variation can 
occur. The interpolation procedure for a general variable, 
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♦ - E N <P (4.1) 

i-1 1 1 

where i is an index which varies with the node number, N 
is the shape function for the i th node and is the value 

^ u 

of the variable at the i node. The Lagrangian shape 


functions for the quadratic, isoparametric element on the 


parent shape are [44] 

N x - - {)(T? 2 - i?) (4.2a) 
N 2 - |(€ 2 + i)(V 2 - V) (4.2b) 
N 3 - j(€ 2 + «)(?? 2 + V) (4.2c) 
N 4 - |(€ 2 - €)(^ 2 + V) (4. 2d) 
n 5 - |(1 - 4 2 m 2 - V) (4.2e) 
N 6 - + 0(1 - V 2 ) (4.2f) 
N? - j(l - t 2 )(V Z + V) ( 4 . 2g) 
W 8 - |(€ 2 - Od " V 2 ) (4.2h) 
N 9 - (1 - * 2 )(1 - i? 2 ). (4.21) 


The Lagrangian shape function, tT, is a continuous func- 
tion having a value of unity at node i and a value of zero 
at node j*i. In the numerical implementation, unknown 
variable distributions will be replaced by an approxima- 
tion with the form of equation 4.1. By this approximation, 
a finite number of parameters may then be solved. 
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4 . 2 Equation Discretization 


Having defined the element, the next stage of numeri- 
cal implementation is to discretize the integral equa- 
tions. The discretization concept in boundary element 
methods is similar to that in finite element methods. The 
concept is that by subdividing the domain of the problem 
into smaller subdomains, or elements, the integrals for 
the complete domain may be approximated on a piecewise 
basis. Furthermore, by enforcing certain boundary condi- 
tions between the subdomains, an approximate solution over 
the entire problem domain may be determined. 


The matrix equations are assembled by satisfying 
boundary conditions at an appropriate number of discrete 
points. The specific points used in this development are 
the nodes of the model. If the given boundary condition 
is a pressure boundary condition, the pressure at that 
point in terms of the boundary source distribution, a, is, 
by equation 3.13, 


* K 

“ /°p U) dB + t 477*^^ 
B 1 


(4.3) 


where p^ is the pressure at the I*'* 1 boundary condition 
* 

point and p^ is the fundamental pressure solution from 
the k “ applied source to the l^ - boundary condition 


point. The unknown boundary source distribution, o, on 
each element is approximated by 
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I a, N 
i-1 


l l 


(4.4) 


utilizing equation 4.1. Furthermore, the boundary integral 
term will be evaluated on a piecewise basis (i.e., on an 
elemental basis). Hence, equation 4.3 can be rewritten as 

N M * K t 

Pi z z s Z OiN.p (r,) dfl + L 4Jry p (4.5) 

1 n-1 B i-1 1 1 1 n k-1 K K1 

n 

where N is the number of boundary elements, M is the 
number of nodes per element, K is the number of applied 
sources, r^ is the variable distance from the 1^ boundary 
condition point to the n th element and B n represents the 
surface of the n^* 1 element. Since the boundary source 
values, cr , are constants, the integral in equation 4.5 
may be rewritten as 

N M * 

* J ( E o Ni)p (r,) dB - 

n-1 B i-1 1 1 ± n 

n 

M N t 

£ a i £ S N iP (r,) dB - Z a a (4.6) 

i-1 1 n-1 B x 1 

n 

The numerical evaluation of the a^ term will be discussed 
in sections 4.3 and 4.4. Finally, equation 4.3 can be more 
simplistically rewritten in approximate form as 


Pl : Z a i a ll + k f 1 4 I ' t k p kl 


(4.7) 


A similar approach can be used for the velocity and 
impedance boundary conditions. With these two types of 


boundary conditions, an additional term due to the singu- 
larity integration is present. To show the consequence of 
the 'free terms', the equation discretization for the 
velocity boundary condition will be presented. In equation 
3.14, the velocity in terms of the boundary source distri- 
bution is given as 


r * Aitco + /ay (r,) dB + I 417*. v. , 
1 B 1 k-1 K K_l 


* 

ki’ 


(4.8) 


Evaluating the integral on a piecewise basis yields 


N 


M 


v^ Z 4 ire 


U°i * n»l l i 'i a ‘ N ‘ V (I 1 ) dB " 
n 


+ E 4^ v v , 
k«l 1 KA 


(4.9) 

where the variables maintain the same meanings as in equa- 


tion 4.5. The 'free term', c^, is zero when i«*l. When 
i«l, c is -j for a point on a smooth boundary and 
for a point on a corner where ft is the solid interior 
angle at the corner. The constant source distributions can 
be factored out of the integral. 

N M * 

4ffc °i ♦ L J ( t o n )v (r ) dB - 

n-1 B l-l 
n 

M N * 

E a (4*c + L J N t v (r,) dB_ = 
i-1 1 ±± n-1 B n 1 ^ n 

E o i (47Tc il + a ijL ) (4.10) 

The final approximation equation for a velocity boundary 
condition is 
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W>£2t^; '■ . V 


V 1 : £ a i( 4rr Cil + a u ) + ^"Vkl < 4 * n > 
Likewise, for impedance boundary conditions, the integral 
equation is 


/ ° 1 (P (*]_) - ZjV ( r x ) ) dB n - 4wz 1 c ll o i 
B n 

K 

+ /l W * (P «- " Z 1 V *1> 1 ° 


(4.12) 


on a piecewise basis or 


K 


Z V 4 " 2 ^!]. + a il> + J^'VPki - Z l V kl> : 0 ^ 4 - 13 > 
for the approximation equation. 


For all boundary conditions, a general matrix equa- 
tion can be written as 


[ A + A' ] 


(4.14) 


where either equation 4.7, 4.10 or 4.13 are used as each 
row of the matrix. The evaluation of the matrix terms for 
each boundary condition type are: 

* 

1. For Pressure Boundary Conditions 


*11 * i N lP < r l> dB n 
n 

a' - 0 
ll 

K 

d k " 4Tr *k p kl 
k*l 


(4.15a) 

(4.15b) 

(4.15c) 


a l * Pi 


(4 . 15d) 


. **» >»* ♦ .ye. 


- . v. . 


2. For Velocity Boundary Conditions 

Si ' l Njv'cri) dB n 
B n 

a ’il " 47rc il 
K 

d v - E 47TF. v ul 
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(4.16a) 

(4.16b) 
(4.16c) 
(4. 16d) 


3. For Impedance Boundary Conditions 

a u “ / ( P*(r 1 ) - z^*^) ) dB n (4.17a) 

B n 

a' iX - -4^2 ^c ^ x (4.17b) 

K * 

d k * k f 1 4 ’ ,, ’k < Pkl - z l v kl > (4-170 

* 0 (4.17d) 

In the D matrix terms, the fundamental solutions are from 
the applied source to the 1 node. The theory can be 

implemented into a computer format using equations 4.14 
through 4.17. 


4.3 Primary Numer ical_ Integration Procedure 

To numerically evaluate the integrals in equations 
4.15a, 4.16a ar.d 4.17a, a numerical integration procedure 
must be utilized. A numerical integration procedure using 
a Gauss -Legendre quadrature was selected. The numerical 
evaluation of the integral equations is performed on a 
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piecewise basis by mapping the integral in the global ele- 
ment coordinate system into the local element coordinate 
system. The domain of the local element is defined between 
the 4 and V coordinates of -1 and 1. The limits are con- 
sistent with those of a Gauss-Legendre quadrature as the 
quadrature formula is [45] 

1 NG 

S g(0 d* : E g(€ i )w i (4.18) 

-1 i-1 1 1 

for a one-dimensional integration of a general function, 
g(£). The values 4;^ and w^ are the integration point and 
weighting factor, respectively, and NG is the total number 
of integration points used in the integration. A boundary 
integration is essentially an area integration (when con- 
sidering three-dimensional problems) which may be 
represented by 

1 1 

i q(t,v) dB - S { g (t,v) d«d7i . (4.19) 

B n -1-1 

For a two-dimensional integration, the quadrature formula 
is [46] 


J g((,V) dB n - E 
B n m-1 


NG , 

L 

n-1 


g(« 


n 


v 

m 



n i i 

Equation 4.20 was used to numerically evaluate 

dary integrals everywhere except where a 

exists in the fundamental solutions as will be 


w. 


m 


(4.20) 


the boun- 
singular ity 
dicussed in 


section 4.4. 


The location of the integration points and the values 
of the weighting factors must be determined to finalize 
the local element integration procedure. A 4x4 Gaussian 


quadrature having coordinates 

<1 - - -0.861136311594053 (4.21a) 
* 2 - V 2 - -0.339981043584856 (4.21b) 
( 3 - V 3 - 0.339981043584856 (4.21c) 
*4 - t ? 4 - 0.861136311594053 (4.21d) 


was selected. The mesh is shown in figure 4.3. The weight- 
ing factors corresponding to the integration points are 


w x - 0.347854845137454 (4.22a) 

w 2 - 0.652145154862546 (4.22b) 

w 3 - 0.652145154862546 (4.22c) 

w 4 - 0.347854845137454. (4.22d) 


The abovement ioned integration quadrature is exact 
for a seventh-order or lower polynomial, g({)- The func- 
tions being integrated in the boundary integrals are the 

* a 

fundamental solutions, p and v times the shape functions 
, tT, which are quadratic polynomials. The fundamental 
pressure solution is 

* e"3 kr 

P ( £ / x ) - -jwp “ (4.23) 


The series expansion of the fundamental solution is 


Figure 4.3 - Integration Point Mesh 
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(4.24) 


Hence, from equation 4.24, some of the limitations of this 
integration can be illustrated. The function will be accu- 
rately integrated only if the first few terms of the 
series approximation dominate. The ^ term poses a problem 
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! 

near a singularity point and will be addressed in the next 

i 

section. 


The preceding discussion applies equally as well for 
the fundamental velocity solution but will not be 
presented here. 


The final point to be made for the integration pro- 
cedures involves transformation from the global domain. 
Since the integration is for the global element, a 
transformation between the domains is necessary as follows 
(47]. 


J*G(x,y,z) dA - Jg(f ,7?) j j(f ,tj) j dfdT) (4.25) 
In equation 4.25, the function G(x,y,z) is a general func- 
tion written in terms of the global element coordinates, 
x, y and z. The function g(f,7j) represents the transforma- 
tion of G(x,y,z) into the local coordinates, i and v . The 
quantity |J(l,i?)| is the determinant of the Jacobian 
matrix For the application at hand, a three- 
dimensional global element is mapped into a two- 
dimensional local element. The determinant of the Jacobian 
for this transformation is [48] 


where 


1 J ( « , i? ) i - ( 


2 

a l 


2 

a 2 


+ a^ 


1 

2 


(4.26) 


- , iv a z 

1 ' 0f 07? 07) 0f ’ 


E 


(4.27a) 






( 


a. 


i* 

a i an ai an 

, is ix _ i!x , 
1 a* an a* an ' 


) 


The partial derivatives in equation 4.27 
with interpolations such as 


lit 

a* 


9 

£ 

i-1 


ae 

3 

i a* 


as an example. 
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(4.27b) 

(4.27c) 

are determined 


(4.28) 


The procedure described above can be applied to any 
element which does not contain a singularity point. When 
the value of r in the fundamental solution approaches 
zero, special considerations are necessary for the 
integration. These considerations are explained in the 
following section. 


4.4 Singularity Numerical Integration Procedure 


A singularity occurs in the fundamental solutions 
whenever the particular boundary condition under con- 
sideration is located on a node of an element over which 
the boundary integral is being evaluated. The variable r 
in the fundamental solutions is the distance between the 
boundary condition point and the integration point. The 
primary Gaussian quadrature is unable to adequately evalu- 
ate the integrals whenever this singularity occurs. Hence, 
a new quadrature had to be developed to handle the singu- 
larity situation. 


In some cases, the singularity can be removed from 
the integration domain. As detailed in the previous sec- 
tion, the integration domain is a quadrilateral region in 
a rectangular coordinate system. However, if the region 
were considered in polar coordinates, the differential 
area is 

dA - rdrdS . (4.29) 

The radius appearing in equation 4.29 may be used to can- 
cel the ^ factor in the fundamental solution. Consider, 
for example, the fundamental pressure solution. 

-jkr . , 

-jup / — dB = -}wp J e~ 3 r drde (4.30) 

B r 3 

n n 

A similar result occurs for the fundamental velocity solu- 
tion . 


I 

B. 


S-v r( 3k )e-^ r 


drde 


(4.31) 


n 


In equation 4.31, it is apparent that a ~ singularity 
still exists even after transforming to a polar coordinate 
system. However, the dot product of the vector normal to 
the surface and the gradient of the radius is zero for a 
flat place element and tends to zero as the radius of cur- 
vature approaches infinity for curvilinear elements. For 
example, in figure 4.4, the fundamental solution is 
evaluated near the upper right corner. The n'^r will tend 
to zero. Thus, the integral in equation 4.31 is integrable 
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aa a Cauchy principle value. (Hence, the 'free term' men 
tioned in chapter 3.) 



Figure 4.4 - Vectors For A Flat Element 


There are three singularity integration cases which 
arise. Those cases are: 

1. Boundary condition at a corner node 


2. Boundary condition at a midside node 
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. . r;-: 

3. Boundary condition at the center node 
One polar domain transformation cam be developed to accom- 
modate all three cases. The transformation for a quadrant 
will apply for a corner node singularity. In addition, the 
remaining two cases cam be considered as multiple qua- 
drants and utilize the transformation for the corner node 
singularity. The transformation can be applied to both 
quadrants neighboring a midside node singularity and all 
four quadrants neighboring the center node singularity. 
Figure 4.5 shows the three cases. The shaded areas 
represent the quadrants transformed to a polar domain. 
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In general, a Gauee-Legendre quadrature similar to 
that explained in the previous section can be used for the 
polar domain. The difference lies in the integration vari- 
ables, in this case r and 9. 


The integration points in the polar coordinate system 
can be transformed from the rectangular coordinate system. 
The values of the <; and V coordinates range from -1 to 1. 
In the circular sector, the radial coordinate ranges from 
0 to 1 and the angular coordinate from 0 to Mapping 
functions can be established which allow the polar coordi- 
nates to be obtained from the rectangular coordinates. 
These mapping functions are 

r - -|(« + 1) (4.32a) 
for the radial coordinate and 


0 « f(7? + 1) (4.32b) 
for the angular coordinate. With equations 4.32, the polar 
integration points become a function of the rectangular 
integration points. 

Determination of the polar weighting factors requires 
more work than for the determination of the polar integra- 
tion points. Weighting factors are evaluated by the equa- 
tion 


b n 
“ / n 
aill 
3*1 




d* 


w . 

i 


(4.33) 
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for a Gauss-Legendre quadrature. To illustrate the deter- 
mination of the polar weighting factors, consider a two 
point quadrature with points and ^ * n the rectangular 
system and points r^ and in the polar system. The equa- 
tion for the first weighting factor in the rectangular 
system is 


w, 


- s 


1 € - i 


-1 *1 " *2 


d£ . 


(4.34) 


For the first weighting factor in the polar system, the 
eauation is 


1 r - r, 


w.. " / 

1 0 *1 


dr 


(4.35) 


Using equation 4.32 with equation 4.35, the weighting fac- 
tor may be expressed as 


w_ 


1 + 1) - * 1} 

-i |(e x + i) - + i) 


(§Sf) 


or 


w 


- if 


+ i) - u, + i) 


(«i + 1) ~ («, + 1) 


d* 


or 


w 


- V 


1 i - €. 


2-1 <1 - *2 


d£ 


Equation 4.38 may now be written as 


(4.36) 


(4.37) 


(4.38) 


w 


- k 

r i 2 «i 


(4.39) 


by substituting in equation 4.34. The 


resulting 


relationship is independent of the number of integration 
points but dependent solely on the multiplicative factor 
in the mapping functions (i.e. r equations 4.32). Conse- 
quently, the weighting factors for the polar domain are 

w r - ^w { (4.40a) 
for the radial weighting factor and 

w 0 - ^w v (4.40b) 
for the angular weighting factor. The weighting factors 
and the integration points in the rectangular domain may 
be found in any reference containing Gaussian quadratures 
C 49 ] . 


In the circular domain, a 4x16 mesh of integration 
points was selected; 4 radial coordinates and 16 angular 
coordinates. The integration points and weighting factors 
for the four point scheme are shown in equations 4.19 and 
4.20. Por the 1£ point scheme, the integration points and 
weighting functions are 


*2 " 


*3 " 

*4 ” 


*6 * 


-0.989400934991649 

-0.944575023073232 

-0.865631202387831 

-0.755404408355003 

-0.617876244402643 

-0.458016777657227 


(4.41a) 
(4.41b) 
(4.41c) 
( 4 . 4 Id ) 
(4. 41e) 
( 4 . 41f ) 
(4.41g) 


-0.281603550779258 
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mfetmd 




’ 8 ' 

-0.095012509837637 

(4 . 41h) 

*9 

a 

0.095012509837637 

(4.411) 

’lO 

m 

0.281603550779258 

(4.413 ) 

’ll 

m 

0.458016777657227 

(4.41k) 

’l2 

m 

0.617876244402643 

(4.411) 

’13 

m 

0.755404408355003 

(4.41m) 

’14 

m 

0.865631202387831 

(4 . 41n) 

’15 

m 

0.944575023073232 

( 4 . 41o ) 

’16 

m 

0.989400934991649 

(4 . 41p) 


W 1 

m 

0,027152459411754 

(4.42a) 

W 2 

m 

0.062253523938647 

(4.42b) 

W 3 

m 

0.095158511682492 

(4.42c) 

W 4 

m 

0.124628971255533 

( 4 . 42d ) 

W S 

m 

0. 149595988816576 

(4. 42e) 

W 6 

m 

0.169156519395002 

( 4 . 42f ) 

W 7 

m 

0.182603415044923 

(4 . 42g) 

W 8 

m 

0.189450610455068 

(4 . 42h) 

W g 

m 

0.189450610455068 

(4.421) 

W 10 

- 

0.182603415044923 

(4.42j) 

W 11 

- 

0. 169156519395002 

(4.42k) 

W 12 

- 

0.149595988816576 

(4.421) 

W 13 

m 

0.124628971255533 

( 4 . 42m) 

W 14 

- 

0.095158511682492 

(4 . 42n) 

W 15 

- 

0.062253523938647 

( 4 . 42o ) 

W 16 

- 

0.027152459411754. 

( 4 . 42p ) 


•i. _-«.\ 
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V 


The mesh for the circular quadrant is illustrated in fig- 
ure 4.6. 


Examination of figure 4.6 reveals a gap between the 
circular domain and the original rectangular domain. As 
shown, the boundary integration for the circular domain 
would exclude that portion of the boundary outside the 
circular domain yet still within the original rectangular 
domain. This omission was corrected with a second domain 
transformation from the circular domain to a rectangular 
domain. The values of the t' and V' coordinates range 
from -1 to 0 since a quadrant of the element is being con- 
sidered. Note that primes are used to differentiate the 
rectangular coordinates in the second transformed domain 
from the original rectangular coordinates. The mapping 
functions for the integration points are 


£' - r - 1 

for 

•<! 

(4.43a) 

£' - rcote - 1 

for 

0 >! 

(4.43b) 


and 

i?' - rtane - 1 for 9 < ^ (4.43c) 

V 1 - r - 1 for 9 > ^ (4.43d) 

For the weighting factors, the mapping functions are 


and 


w, , - w 
£ ' r 


(4.44) 


w. 


V - 


for » < ^ 


cos® 


(4.45a) 




M " \S ' 
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w 


0 


V • ,in9 fot 8 > 4 <4. 45b) 
These mapping functions were determined in the same manner 

as the mapping functions for the first domain transforma- 
tion (i.e., rectangular to polar). Figure 4.7 shows the 
mesh for the transformed circular domain. 



The final term required for the numerical integration 
ie the determinant of the Jacobian. For a normal integra- 
tion procedure, the determinant of the Jacobian is com- 
puted using equations 4.26 and 4.27. The determinant of 
the Jacobian can also be viewed as a ratio of the global 
area to that of the local area. 


dxdy 

d£di? 


(4.46) 


In polar coordinates, equation 4.46 is rewritten as 


i w i 


r dr d® 

-3-9.JL 

r^dr^d®^ 


(4.47) 


where the subscripts indicate the global and local domain 
variables. For the singularity integration, the integral 
is of the form 


h ( V e g )dr g d0 g " / F Ui* V I J l r dr l d ®l (4.48) 
where i J S r is the determinant of the Jacobian required for 
the singularity integration procedure. From equation 4.48, 
the required determinant of the Jacobian will be 


drdfl 
, j, _ _9 g 
,J, r dr-jd^ 

and can be calculated from the normal determinant 
Jacobian by 


(4.49) 
of the 



9 

If the mapping between the global and local domains is 
simply a scaling of the same angular dimension sector, the 


mapoing equations can be written as 


9 

for the radial component and 


ar 


(4.51a) 


« g " 9 i (4.51b) 
for the angular coordinate since it is assumed that the 
element is only varying in size and not in shape. Thus, 
the differentials are expressed as 


dr - adr, 
9 1 


(4.52a) 


and 


d®g - d0 x . (4.52b) 

For this case, the normal determinant of the Jacobian is 


IJ{ - 


r_dr_d0 
9 9 9 

r 


(ar ) (adr )d6 
1 11 

r l dr l dS l 


(4.53) 


and the required determinant of the Jacobian is 


dr _d0 

—3 2 

dr l d0 l 


(adr^d^ 

dr l de i 


- a 


(4.54) 


Therefore, the required determine - t of the Jacobian may be 
calculated from the normal determinant by 




nJ737 


(4.55) 


The relationship in equation 4.55 is an approximation if 
the shape of the global element is distorted from that of 
the local element. Nevertheless, this relationship is 
utilized for the singularity integration. Hence, the qua- 
drature formula becomes 
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mat -V^'v ■•%.•-. 


Sf (x,y)dxdy - /F (r ^ 9^ , j ; dr jd01 . (4.56) 

The numerical integration procedures detailed above 
car. be utilized to evaluate the discretized integral equa- 
tions to establish the system matrices for a boundary ele- 
ment model. Although this numerical integration procedure 
is successful, it is not the only possible procedure to 
apply. An alternative integration proc«dure has been sug- 
gested by Zimmer le and Bernhard [50] but was not attempted 
in this research. 

4.5 Solution Procedure 

Equation 4.14 is used to solve for the assumed source 
densities once the numerical procedures described above 
have been used to evaluate and assemble the system 
matrices. To solve for the source densities, a full 
matrix, complex equation solver is required. The matrix 
equation is neither banded nor symmetric so a robust equa- 
tion solver is required. The subroutine used, the IMSL 
subroutine LEQT1C, is a linear equation solver for equa- 
tions having complex quantities. It uses a full matrix 
storage mode. This storage mode is required since tlv 
matrices for boundary element methods, unlike finite ele- 
ment methods, are fully populated matrices. The IMSL sub- 
routines are readily accessible as a co>nmercial code which 
was the primary reason for the application in this 
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research. It is possible that a more efficient equation 
solving routine is available or may be generated for use 
in solving the boundary integral equations. 

The goal for this research was to predict the sound 
fields in acoustical cavities. To accomplish this, the 
source densities, after quantization, were used to solve 
for che acoustic pressures at field points interior to the 
cavity domain. Equations 4.14 and 4.15 were again util- 
ized. The unknown variables, the pressures at the field 
points, can now be solved utilizing the source density 
distribution found earlier. The only major difference in 
the formulation is t\ - ; :he distance r in the fundamental 
pressure solution is the distance from the source density 
to the field point at which the pressure is desired. Using 
this procedure, the pressures at the interior field points 
may be evaluated. 

The procedures required for numerical implementation 
of the boundary integral equations have been discussed. To 
illustrate the techniques, a number of examples were con- 
sidered. The numerical results of these examples are 


presented 


in the following chapter. 


CHAPTER 5 


CASE STUDIES 


A number of case studies were undertaken to illus- 
trate the theory and numerical implementation discussed in 
chapters 3 and 4. Two models were considered in these stu- 
dies: a spherical cavity model and a rectangular cavity 
model. This chapter presents the numerical and, where 
available, the analytical results of the case studies. 

5.1 Spherical Cavity Model 

A spherical cavity was chosen as the first case study 
for a number of reasons. A sphere with a velocity distri- 
bution which is uniform with respect to angle can be 
treated as a one-dimensional problem as will be demon- 
strated in the derivation of the analytical solution. 
Consequently, the analytical solution for the spherical 
cavity can be easily determined. Another advantage for the 
sphere is the absence of corners and edges. Corners and 
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create difficulties with boundary element methods as 
>e shown for the rectangular cavity. 


’he analytical solution for the spherical cavity is 
;d from the wave equation for spherically symmetric 
ire fields [51]. 


-,2 . 2 
J . i Le 

-2 ' c 2 at 2 


a r ‘ 


(5.1) 

tng the time dependence from equation 5.1 yields the 
Lon 


a p 2 8p 

7 + - + 

a t 2 1 dr 


k 2 p 


(5.2) 


represents the time-independent wave equation in a 
ically symmetric geometry without applied sources, 
equation is a special case of the spherical Bessel's 
Lon 


r 2 ^ + 2r|^ + k 2 (r 2 - n 2 )p - 0 (5.3) 

2 

ed by r and with n-0. Thus, the solution to equation 
l 11 involve spherical Bessel functions and will be of 
5rm 


p(r ) - Aj 0 (kr) + By o (kr) (5.4) 

the constants A and B are determined by the boundary 
r .ion8. Alternatively, a second general solution to 
oherical Bessel's equation of equal validity is 

„ -jkr _ jkr 

P(r ) - a ' + b ^ — (5.5) 


where a and b are determined by the boundary conditions. 


For the first verification of the numerical code, the 
specific solution of 


p(r ) - 100j Q (kr) (5.6) 
was selected. The spherical model has a radius of meters. 
Also, am arbitrary wave number of 8 was selected. This 
corresponds to a frequency of 437 Hz in air with a speed 
of sound of 343 meters/second. Figure 5.1 shows the 
analytical solution of the pressure distribution within 
the the cavity. The sphere center corresponds to a radial 
position of zero while the sphere boundary corresponds to 
a radial position of 


Magnitude Phase 
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Figure 5.1 


Analytical Spherical Pressure Distribution 
For k-8 
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Appropriate boundary conditions of various types were 
derived to verify the solution in the cavity interior. The 
solution of equation 5.6 at the boundary is p(^) - -18.92 
with the selected wave number. For the first example, 
this pressure was applied at the outer boundary of the 
sphere. The numerical results are shown with the analyti- 
cal results in figure 5.2. The numerical solution is 
represented by the solid line marked with Y's at the 
discrete solution points. The dashed line represents the 
analytical solution. 

A sphere with uniform velocity boundary conditions 
was the second numerical example. The velocity boundary 
conditions were determined by 

(5.7) 

which relates the velocity of a spherical wave to the 
pressure [51]. For the numerical example, the velocity 
boundary condition was v - -0.0456 + j0.0114. Figure 5.3 
8 hows the numerical and analytical results. 
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Figure 5.3 - Spherical Pressure Distribution For Velocity 
Boundary Conditions At k-8 
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The analytical and the numerical results in both fig- 
ures 5.2 and 5.3 are in excellent agreement. In fact, the 
correlation is so good for figure 5.2 that the two results 
are nearly coincident. One note must be made regarding the 
phase plot. When examining the phase, it must be under- 
stood that a phase of 180° is equivalent to a phase of 
-180°. Hence, the apparent jumps in phase in the region 
from r » 0.4 to r - 0.5 do not actually represent phase 
differences as it may appear at first glance. For the 
velocity boundary condition model represented in figure 
5.3, a slight amount of difference can be detected between 
the analytical and numerical solutions. This slight error 
is due to the interpolation of the radial distances 
between the boundary condition points and the integration 
points. As previously noted, the shape functions are qua- 
dratic functions. However, the radial distance is deter- 
mined by 


r " -J(* 2 " x l )2 * (y 2 ” y l )2 + (z 2 " 2 1 )2 (5 - 8) 

# 

which is not a purely quadratic function due to the radi- 
cal involved in its determination. Consequently, some 
interpolation error can be expected and, in fact, is evi- 
dent. The interpolation error is more prevalent when the 
velocity boundary conditions are used because of the ~ 

term in the fundamental velocity solution as opposed to 


the term in the fundamental pressure solution. Regard- 
less, for cases where derived boundary conditions are 
specified, the prediction of the solution in the domain is 
very good. 


The next case study was that of a pulsating spherical 
cavity. The analytical solution for this problem can be 
represented by equation 5.5 where the boundary conditions 
for this problem are 


_ 

dr 


0 


&'t the center of the sphere and 


(5.9a) 


- - jwpv(r Q ) (5.9b) 
at the surface of the sphere. Boundary condition 5.9a 
respresents the condition that the pressure is finite at 
the sphere center. Equation 5.9b is derived from Euler's 
equation. Using equation 5.5 with the boundary conditions 
of equations 5.9, the analytical solution for the pressure 
is 


pu) - y 


P c * r o v(r o> 


(sin(kr Q ) - kr Q cos (kr Q ) ) 


sin(kr ) 
r 


(5.10) 


For the model in question, the radius at the sphere sur 


face is r - 

o z 


Resonant frequencies can be determined for the spher- 
ical cavity using equation 5.10. Resonance occurs whenever 
the impedance ^ becomes infinite. From equation 5.10, 


resonance for the spherical cavity occurs whenever 

sin(kr ) - kr cos(kr ) =* 0 . (5.11) 

o o c 

The transcendental equation for the natural frequencies is 

tan(kr Q ) - kr Q (5.12) 

for the spherical cavity. Thus, the first three natural 

frequencies occur at kr - 0, kr - 4.493 and kr - 7.725. 

o o o 

The analytical solution of equation 5.10 at r - r Q is 
shown in figure 5.4 as a function of the frequency. The 
natural frequencies are evident in figure 5.4 as peaks in 
the pressure magnitude. In addition, figure 5.4 shows zero 


pressures 

and 

180° phase 

shif ts 

occur ing at 

kr Q - tt and 

Kr o 

- 2^. The zero 

pressures 

are expected 


since p ■ 0 whenever sin(kr Q ) - 0 in equation 5.10. The 
velocity at the boundary of the sphere was chosen to be 1 
to generate the analytical solution in figure 5.4. 

To numerically evaluate the response of the spherical 
cavity, a sweep of frequencies from 0 to 873 Hz was used 
to generate results comparable to those in figure 5.4. The 
frequency range corresponds to kr Q values from 0 to 8. A 
uniform velocity boundary condition of 1 was applied to 
the spherical cavity model. Both the analytical and numer- 
ical results are shown in figure 5.5. The solid line 
marked with Y'e represents the numerical results and the 
dashed line is the analytical solution. 
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To further investigate the spherical cavity response, 
pressure distributions within the cavity were numerically 
determined at three frequencies. The three frequencies are 
kr Q • 2, 5.4 and 7.1 corrseponding to 1372, 3704 and 4871 
Hz, respectively if r Q ■ meters and c - 343 m/s. The 
spherical cavity pressure distributions for each of these 
frequencies are shown in figures 5.6, 5.7 and 5.8. All 
three distributions exhibit the expected spherical 
Bessel's function-type behavior. The particular frequen- 


cies 

were 

selected 

to show the 

behaviors 

occur ing 

before 

and 

after 

the 

zero pressure frequencies. 

i.e. , 

kr - 

v and 

2n. The 

analytical 

solutions 

shown in 

these 


figures are from equation 5.10. As in figures 5.2 and 5.3, 
excellent correlation is shown in figures 5.6 through 5.8. 
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Figure 5.6 - Spherical Cavity Pressure Distribution For 
Velocity Boundary Conditions At kr-2 
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Figure 5.7 - Spherical Cavity Pressure Distribution For 
Velocity Boundary Conditions At kr-5.4 




Magnitude 


4000.0 

3000.0 

2000 .0 
1000.0 
0.0000 

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .45C 

Radial Position 


Figure 5.8 - Spherical Cavity Pressure Distribution r r 
Velocity Boundary Conditions At kr=7.1 








89 




Although the results presented thus far have only 
been for pressure and velocity boundary conditions, 
impedance boundary conditions have also been implemented. 
Models having only impedance boundary conditions must con- 
tain an applied source within the cavity. Otherwise, the a 
vector in equation 4.14 is a zero vector (see equation 
4.17d) and the resulting source distribution is trivial. 
In other words, there is no excitation of any sort in the 
problem. To verify the impedance boundary condition imple- 
mentation, a spherical cavity with purely impedance boun- 
dary conditions and an applied source of strength 1 
located at the sphere center was utilized. The impedance 
boundary conditions were established such that they would 
simulate the impedance experienced by a point source radi- 
ating into an infinite free space. The impedance relation- 
ship can be derived from equation 5.6 as 


z - 


Pc 


d - Bf> 


(5.13) 


With a wave number of k * 0.583 (200 Hz) and r * t at • the 
boundary, the resulting impedance is z - 32.5 + jlil.5. 
Figure 5.9 is the pressure distribution obtained numeri- 
cally for the described model. 
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Foe comparative purposes, similar distributions for 
models with pressure and velocity boundary conditions were 
obtained. The pressure solution and, consequently, the 
boundary conditions may be determined in a number of ways. 
In this case, the pressure boundary conditions were deter - 
mined from equation 3.13 and 3.15. 


— * — 


* — — 


P(x) - Jo(Op (C,x) + t 4wt.p U ai ,x) (5.14) 
B i-1 1 1 

The pressure at any point is due entirely to the applied 
source if the boundary simulates an infinite domain prob- 
lem. Hence, the source density distribution is zero at 
the m^del boundary (i.e., o(C) * 0). Equation 5.14 becomes 


e -j*r 

p(r ) - -j<i>p4^ 1 — - — (5.15) 
after the fundamental pressure solution is inserted. With 
an applied source strength of 1, the equation to determine 
the pressure boundary condition is 


p(r) - -jwp4fr — - — . (5.16) 
For ths model with k - 0.583 and r - j at the boundary, 
the pressure boundary condition is p - -139.07 - j463.51. 
The numerical results obtained are shown in figure 5.10. 
Equivalent velocity boundary conditions found using equa- 
tion 5.6 are v - -4.167 + j0.G33. Figure 5.11 shows the 
pressure distribution numerically obtained using velocity 
boundary conditions. 
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To evaluate the validity of the impedance boundary 
condition capability of the program, figure 5.9 may be 
compared with figures 5.10 and 5.11. Both pressure and 
velocity boundary conditions have previously been shown to 
produce excellent results in the spherical cavity model. 
Figures 5.9, 5.10 and 5.11 are identical pressure distri- 
butions. Regardless of the type of boundary condition 
being applied, the numerical results for this model are 
unchanged. 

The final study using the spherical cavity was tc 
numerically generate a forced response prediction for the 
case of a point source in a rigid-walled cavity. Velocity 
boundary conditions of zero were applied to the sphere. A 
source having a strength of 1 was placed at the center of 
the cavity to excite the model. The swept frequency band 
ranged from 0 to 873 Hz (or kr from 0 to 8) as was the 
case for the previous response in figure 5.5. The response 
for the rigid-walled spherical cavity is shown in figure 
S.12. In addition, a pressure distribution within the cav- 
ity for a frequency of 3087 Hz (kr - 4.5), near the first 
resonant frequency, is shown in figure 5.13. 


(X10 > Phase 
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The rigid-walled spherical cavity response compares 
well with the response in figure 5.5 in terms of the 
resonant frequency identification. Both responses identify 
peaks occuring at approximately 3080 Hz and 5300 Hz (kr's 
of 4.493 and 7.725). For the rigid-walled cavity, the 
pressure magnitudes are much higher due to the effect of 
the applied source. The phase relationships differ between 
the two responses for the same reason. The cavity pressure 
distribution on figure 5.13 exhibits the expected Bessel’s 
function-type behavior. However, the effect of the applied 
source can be seen near the radial position of zero. Since 
the pressure at the applied source is infinite, the near- 
field effect close to a radius of zero is expected. Other- 
wise, the pressure distribution is dominated by the 
response of the cavity near the resonant frequency. This 
behavior could be expected to be much different if the 
cavity was being excited at a frequency away from a 
resonant frequency. Note that the low frequency regime in 
figure 5.12 is characteristic of the low frequency regime 
in figure 5.5. Once again, the low frequency results are 
affected by the formulation and do not approach infinity 
as they should. 

In general, the indirect boundary element method has 
done an excellent job for the variety of spherical cavity 
problems tested. The models for boundary excitation, 
internal point source excitations and with impedance 
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boundary conditions all perforin very well with the one 
exception at low frequency. 

5.2 Spherical Model For Radiation 

Although the research being presented deals primarily 
with cavity acoustics, it is possible to deal with acoust- 
ical radiation problems applying the same formulation. To 
do so, the models need only have negated 'free terms' 
(i.e., the quantity c in equations 4.16b and 4.17b). In 
other words, instead of having the 'free terms' be nega- 
tive (as is done for cavity problems), the 'free terms' 
must be positive. 

The acoustical radiation model for the sphere was 
studied only to determine whether the formulation would 
fail at the interior eigenf requencies for the exterior 
radiation problems. As mentioned in the previous section, 
this type of formulation failure has been demonstrated by 
other researchers. 


Radiation from a pulsating sphere can be modeled by 
applying uniform velocity boundary conditions to the 
sphere and applying the 'free term' negation. The analyt- 
ical solution for the pressure on the surface of a pulsat- 
ing sphere is [52] 


kr 


p(r n ) - PcV 


°((kr Q ) 2 + 1) 


(*r 0 + 7) 


(5.16) 


99 


Equation S.16 was utilized to generate the analytical 
boundary pressure response in figure 5.14 for a velocity 
boundary condition of 1. 

The sphere model was used to numerically determine a 
boundary pressure response plot for a pulsating sphere. A 
uniform velocity boundary condition of 1 was applied to 
the sphere. The 'free terms' were negative of those used 
for the spherical cavity studies. Figure 5.15 contains the 
analytical solution represented by the dashed line and the 
numerical solution represented by the solid line marked 
with Y's at discrete frequency points. 
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Figure 5.14 - Analytical Pulsating Sphere Boundary Pres 
sure Prediction For Radiation 
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Pigure 5.15 vary clearly shows that the formulation 
failures for the radiation problem are present in the for- 
mulation being applied in this research. At the interior 
eigenfrequencies of kr equal to if and Iff, the boundary 
pressures diverge from the analytical solution. Koopmann 
and Benner [29] present results for a pulsating sphere 
which portray the same behavior as shown in figure 5.15. 
Thus, the conclusion to be drawn is that an overdetermina- 
tion procedure should be applied if this formulation were 
to be used for radiation problems. Also, these results 
lend support to the hypothesis that the difficulties 
obtained in the low frequency regime for the cavity prob- 
lems may be due to a formulation failure. Near a frequency 
of zero, rigid-body modes exist which this formulation is 
evidently incapable of determining for cavity problems. 

5.3 Rectangular Cavity Model 

All the results presented heretofore in this chapter 
have been for a spherical model. For boundary element 
methods, corners and edges are known to create difficul- 
ties for boundary element techniques due to the boundary 
discontinuities. These difficulties are described in the 
literature in some detail. Banerjee and Butterfield [12] 
have devoted an entire chapter in their text to the sub- 
ject. Although the problem is well-documented, there does 
not exist one ultimate solution for all cases. 


To investigate the consequences of corners and edges, 
a rectangular cavity with dimensions of 4x5x3 was studied. 
Three different methods of modeling the corners and edges 
were attempted. The first method was to model the corners 
and edges as intersections of adjacent sides. Essen- 
tially, the first method is the intuitive aproach one 
would assume having no knowledge of the inherent difficul- 
ties. The normal vectors for the nodes on the corners and 
edges were given components perpendicular to each inter- 
secting wall. The second method involved leaving gaps 
between the intersecting walls. In essence, the nodes at 
the corners and edges were not physically attached to 
adjacent walls of the box. Hence, compatibility of adja- 
cent node points is not enforced. (Element formulations of 
this type are referred to as ' noncompatible ' elements.) 
The final method involved rounding off the corners and 
edges to avoid the boundary discontinuities. 

To evaluate the three methods of modeling the corners 
and edges, one wall of the cavity was driven while the 
remaining five walls were rigid. A simplistic analytical 
approach for this situation is to consider the problem on 
a one-dimensional basis as a driven/closed cavity. The 
pressure solution for this problem is 

p . A e 5*< L - *> * I e -3 k(L ' x > (5.17) 
where L is the length of the cavity in the direction 
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normal to the driven end, x is the length variable ranging 
from 0 to L and the constants A and B are determined by 
the boundary conditions. The boundary conditions are that 
the velocity at the closed end is zero (V(L) - 0) and the 
velocity at the driven end is equal to the driving velo- 
city (V(0) - V ) . With these boundary conditions and equa- 
tion 5.17, the approximate analytical pressure solution is 

cosk(L - x) 

p - ~J pcV o a in()cL) • < 5 - 18 > 

To quantify equation 5.18, a driving velocity of 1 and a 

dimension of L « 5 (corresponding to the y dimension of 

the cavity) were used. The pressure distribution along 

the y coordinate for wave numbers of k equal to 1, 2 and 3 

(or frequencies of 343, 686 and 1092 Hz with a speed of 

sound of 343 m/s) are shown in figures 5.16, 5.17 and 

5.18, respectively. 
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Pigure 5.16 - Analytical Pressure Distribution For A 
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Boundary element model* incorporating the three 
corner and edge methods were used to numerically evaluate 
a driven/closed cavity. However, the driving velocity was 
not uniform as assumed in the analytical solution. An 
attempt to produce such a situation is impossible in the 
first and third models because the normal velocity at the 
corner nodes must be compatible on both walls. Imposing 
such a boundary condition would cause the adjacent walls 
to be non-rigid in the fii3t and third models. Conse- 
quently, a velocity profile having a magnitude of 1 at the 
center but varying linearly to zero at the edges was 
applied to the driving wall (i.e., the wall at y equal to 
zero). Pressure distributions were calculated for wave 
numbers of k equal to 1, 2 and 3 just as for the analyti- 
cal solution. The numerical results are shown in figures 
5.19, 5.20 and 5.21 for the first model (the true box), 
figures 5.22, 5.23 and 5.24 for the second model (the 
gapped box) and figures 5.25, 5.26 and 5.27 for the third 
model (the rounded box). 
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Figure 5.27 - Pressure Distribution For Third Box Model At 






Two points must be made regarding the numerical 
models. Firstly, the excitation at the driving wall is not 
uniform as in the analytical solution. Secondly, the 
numerical models are three-dimensional. Thus, one can not 
expect purely one-dimensional results, especially when the 
excitation at the driving wall is not purely one- 
dimensional. These points need be highlighted to provide 
explanations for the differences between the analytical 
and numerical solutions. 

With the previous points in mind, some qualitative 
judgements can be made on the three numerical models. The 
analytical solutions in figures 5.16 through 5.18 can be 
used as guidelines to evaluate the three models. For all 
three frequencies, the analytical solutions are 
cosinusoidal in nature. The numerical results for the 
first model appear to exhibit some form of sinusoidal 
behavior for the second and third wave numbers but not as 
clearly defined as would be expected either in magnitude 
or phase. The second model shows more clearly cosinusoidal 
behavior in the pressure distribution. Phase shifts of 
180° are evident at positions of zero pressure magnitude. 
There is some variability in the peak magnitudes within an 
individual pressure distribution. However, it is quite 
possible that the variability is due to cross-coupling 
with modes in the other two orthogonal directions since 
the absolute maximums occur near the center of the cavity 
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(e.g., figure 5.24). Finally, the third model, the 
rounded corner model, shows cosinusoidal behavior at the 
first wave number but not at the second or third wave 
numbers. The pressure magnitudes are not regular peaks nor 
are any clearly defined phase shifts evident. As a result 
of this numerical study, the second method of modeling the 
corners and edges, that of providing gaps, appears to be 
the best of the three methods and is also the one often 
recommended in the literature (e.g., [12]). It should 
also be noted that another solution which should be inves- 
tigated for acoustical problems is the use of ' noncompati- 
ble ' elements. In 'noncompatible' elements, the node 
points lie in the interior of the element rather than on 
the edges. Many of the problems encountered here will not 
be a consideration with 'noncompatible' elements. 

After determining that the gapped edge model provided 

the best results of the three models studied, the gapped 

0 

edge model was utilized to study the internal source exci- 
tation problem of the rectangular cavity. The boundary 
conditions for the response study were rigid-wall boundary 
conditions. Excitation of the model was supplied by a 
source placed in the center of the rectangular cavity. The 
strength of the source was 1. For a rigid-walled rec- 
tangular cavity, the pressure distribution in terms of 
modal response will be of the form [53] 
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p * C Z1 H A^„COS (t )C08(^ )C08(^ ) (5.19) 
1 m n x y z 

where L , L and L are the dimensions of the rectangular 
x y z 

cavity in the x, y and z directions, respectively and A lmn 
is the modal amplitude. The modal amplitude will depend on 
the closeness of the driving frequency to a natural fre- 
quency and also on the coupling of the exciter to the 
mode. Since an applied source is placed at the center of 
the cavity for excitation purposes, the odd modes will not 
be excited. (Odd modes means that the variables 1, m or n 
in equation 5.19 are odd numbers.) The applied source 
creates pressure antinodes in the center of the cavity. 
Hence, only even modes should be excited. 

The numerically determined response plots for the 
rigid-walled rectangular cavity are shown in figures 5.28, 
5.29 and 5.30. The pressures shown in the response plots 
are associated with a point located in the center of a 
wall of the cavity. For example, the x direction response 
is for a point located in the center of a wall having a 
constant x coordinate (i.e., a wall parallel to the yz 
plane). The response plots indicate resonant frequencies 
by the significant peaks in the magnitudes and shifts in 
the phase. Incidentally, the difficulties in the low fre- 
quency regime are apparent again in these response plots 
below a wave number of about 0.25. 








Magni tud 














Pressure distributions within the rectangular cavity 
were generated from the numerical results to identify the 
pressure response shapes occur ing at resonant frequencies. 
Near the resonant frequencies the response should be dom- 
inated by the respective mode shape. Between resonance, it 
wot’ld be expected that the response would be lower and 
would be a superposition of all the near-frequency 
resonant mode shapes. Since the cavity and pressure dis- 
tributions are three-dimensional, the pressure distribu- 
tions were determined for planes perpendicular to one of 
the orthogonal directions and passing through the center 
of the rectangular cavity. Thus, for each response at 
which a distribution is determined, three plots are gen- 
erated; one in each of the three orthogonal planes. Furth- 
ermore, the magnitude at the center of all the magnitude 
plots is infinite and the phase at the center of the phase 
plots is -90° because of the applied source. As a result, 
the response of the center point is truncated to the level 
of the maximum pressure magnitude. 

From the cavity response, the first resonant fre- 
quency occurs at a wave number of 1.34 (460 Hz). The pres- 
sure distribution within the rectangular cavity was gen- 
erated at that wave number. The pressure distribution of 
the variations in the z direction is shown in figure 5.31. 
Likewise, figure 5.32 is the y variations and figure 5.33 
is the x variations. Figure 5.31 indicates a possible 
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cosinusoidal variation in magnitude along the y coordinate 
which is confirmed by the clearly defined 180° phase shift 
along the y coordinate. Relatively constant pressures are 
displayed in figure 5.32. The cosinusoidal variation in 
the y coordinate is evident again in figure 5.33. It is 
readily apparent that the (0,2,0) mode is the resonant 
mode shape at 460 Hz. Analytically, the (0,2,0) mode shape 
occurs at k * or 432 Hz. The conclusion is that the 
boundary element model successfully located the first 
resonant frequency with about, a 6% upward shift in the 
frequency. 







Figure 5.32 - Rectangular Cavity Pressure Distribution At 
k-1.34 (XZ Plane) 
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The second resonance condition occurs at a wave 
number of about 3.65 (566 Hz). The pressure distributions 
within the rectangular cavity are given in figures 5.34 
through 5.36. A cosinusoidal variation along the x direc- 
tion is indicated by the magnitudes and by the 180° phase 
shifts in figures 5.34 and 5.36. Hence, the (2,0,0) mode 
shape appears to be the dominant mode at 566 Hz. The 

2 it 

(2,0,0) mode occurs analytically at k - or 539 Hz. For 
the second resonant condition, the boundary element model 
located the resonance with a 5% upward shift in the fre- 


quency. 
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The third resonant pressure distribution is at a wave 
number of 2.06 (707 Hz). In the response of figures 5.28 
through 5.30, two closely-spaced peaks appear between the 
wave numbers of 2 and 2.25. The first peak is the one at 
707 Hz. Pigures 5.37 through 5.39 are the pressure distri- 
butions of the third resonant condition. In figure 5.37, 
it seems as though a cosinusoidal variation exists along 
both the x and y directions. However, it is not entirely 
clear from the phase relationship that this is indeed the 
case. By examining figure 5.38, it becomes more obvious 
that just such a variation is present in the x direction. 
Similarly, the y direction variation is better viewed in 
figure 5,39. The resulting conclusion is that the (2,2,0) 
mode is the dominate mode at this frequency. Theoreti- 
cally, the (2,2,0) mode occurs at a wave number of 

v (^)^ + (5“")^ * 2.01 (689 Hz). For the third resonance 
in the rectangular cavity response, the boundary element 
model located the resonant frequency with only a 3% upward 


shift. 
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Figure 5.38 - Rectangular Cavity Pressure Distribution At 
k«2.06 (XZ Plane) 
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5.39 - Rectangular Cavity Pressure Distribution At 
k-2.06 ( ZY Plane) 
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Although the pressure distribution will not be 
presented for the resonant peak at the wave number of 2.21 
(758 Hz), some comments can be made regarding it. The next 
expected mode shape should be the (0,0,2) mode which, 
theoretically, is excited at a wavs number of 2.09 (717 
Hz). The numerically determined pressure distribution at 
758 Hz was characteristic of the (0,0,2) mode but also 
retained the nature of the (2,2,0) mode which occured at a 
slightly lower frequency. Because the two resonances are 
so closely related in terms of frequency, both the (2,2,0) 
mode and the (0,0,2) mode are apparent in the distribution 
at 758 Hz. A situation such as this is likely to happen 
whenever the modal density is high in a particular fre- 
quency range cf interest. 

A series of pressure distributions were gensrated for 
a non-resonant frequency of the rigid-walled rectangular 
cavity. The distributions in figures 5.40 through 5.42 are 
for a wave number of 1 (343 Hz) which is a frequency lower 
than the first resonant frequency of 460 Hz. The pressure 
distributions plainly show that at this non-resonant fre- 
quency, the applied source is the sole determinant of the 
pressure field characteristics. This type of distribution 
should be expected as all the resonant response will be 
low at this frequency. 
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5.4 Summary 

The first section of this chapter presented numerical 
results and, analytical results where applicable, of sound 
fields in a spherical cavity. Excellent correlation was 
obtained except in the low frequency regime. A hypothesis 
for the low frequency difficulties was presented and par- 
tially supported through literature references. Using the 
spherical cavity model, all three boundary condition types 
(pressure, velocity and impedance) and the applied source 
capabilities were verified. In addition, the sphere model 
was used to briefly examine the acoustical radiation pred- 
iction capabilities of the program for a pulsating sphere. 

The remainder of the chapter presented results 
obtained for a 4x5x3 rectangular cavity. Three different 
techniques for modeling corners and edges were presented. 
The best of the three modeling techniques was to model the 
corners and edges by leaving slight gaps. A model of the 
cavity having gaps at the corners and edges was utilized 
to develop a cavity response and pressure distributions 
both on and off the resonant frequencies. The resonant 
frequencies were located and identified with good accuracy 
using the boundary element model. 


CHAPTER 6 


CONCLUSIONS 


The objective for this research was to develop a 
numerical tool based on indirect boundary element theory 
to predict the sound field within an acoustical cavity. In 
addition, two additional features were desired. The first 
additional feature was that the three boundary condition 
types; pressure, velocity and impedance, be formulated and 
implemented in the development. Secondly, the implementa- 
tion of acoustic point sources was desired. 

A code to implement these objectives was developed 
and verified. In general, the conclusions which can be 
drawn from this experience are? 

1. The indirect boundary element method is capable of all 
analyses for which direct boundary element methods have 
been demonstrated in the literature. In addition, IBEM 
techniques require only one boundary integral in most 
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cases. This will be more efficient than DBEM's in man/ 
instances . 

2. Forced response in three-dimensional cavities is a 
straightforward application of boundary element 
methods. Corners and edges cause difficulties as they 
do with radiation problems. Further conclusions in this 
regard are included in the next paragraph. 

3. Implementation of internal source capabilities is rea- 
sonably straightforward and potentially quite powerful 
in future applications such as active noise control. 
Boundary element methods are decidedly better then fin- 
ite element methods for representing internal point 
sources since BEM's can accommodate singularities. 

4. The application of impedance boundary conditions was 
successfully demonstrated. It should be pointed out 
that the formulation used in this development is not 
limited to locally reacting impedance. 

One suggestion can be made regarding the problem with 

t 

corners and edges evidenced in the rectangular box model. 
The application of a ' noncompatible ' element may resolve 
these difficulties. A 'noncompatible' element is one in 
which the nodes of the element are not located along the 
edges of the element as for the isoparametric element used 
in this research but located internal to the element. By 
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locating the nodes, or the locations of the boundary con- 
ditions, within the element, the boundary conditions will 
never be situated at a boundary discontinuity. Thus, the 
'free terms' in the velocity and impedance boundary condi- 
tion equations will always maintain a value of and, 
more importantly, the difficulties with the boundary 
discontinuities can be avoided. 

The boundary element procedure developed for this 
research can be applied to noise source identification. 
Although this concept was not specifically pursued in the 
case studies, the identification of acoustical generators 
and optimization of the cavity acoustical characteristics 
can be guided with the application of this boundary ele- 
ment procedure. 

In summary, it appears from this work that the 
indirect boundary element method should be the numerical 
method of choice for a large precentage of investigators 
interested in forced response problems in both interior 
and exterior domains. This is especially true when only a 
limited number of frequencies and response locations are 


desired . 
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